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We have initiated a programme to compute the lower moments of the unpolarised and polarised deep inelastic 
structure functions of the nucleon in the quenched approximation. We review our progress to date. 



1. INTRODUCTION 

Lepton-nucleon deep inelastic scattering (DIS) 
can be described in terms of 4 structure functions 
- Fi, F2, gi, g2- Unpolarised beams lead to Fi, 
F2 which have been studied for many years. Po- 
larised beams are needed for gi and 32 • ffi has on- 
ly recently been measured by the EMC collabora- 
tion [1], while experiments for g2 are planned e.g. 
at HERMES. From the Wilson operator product 
expansion, moments of the structure functions are 
related to certain nucleon matrix elements (see 
e.g. [2]). Thus calculating these matrix elements 
leads to a knowledge of the structure functions. 
For the unpolarised structure functions we have 
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(with even n starting at 2), where the expansion 
is in powers of l/Q^ and 

■•^">|p,s-} =2i;„[p''^ •••p^" -Tr] (2) 
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with 



...£)f"^-Tr. (3) 



Here ip = u or d (quark flavour /, charge Q^^^)- 
The additional Wilson coefficient has been set 
equal to 1 (i.e. leading order in perturbation 
theory). (The normalisation is {p,s\po,So) = 
{2Tif2E^8{p-po)8s,s„ with = -^2.) The 
moments, eq. (1), have a parton model inter- 
pretation, being the powers of the fraction of 
the nucleon momentum carried by the parton: 
Vn = (x"^^). For the polarised structure func- 
tions, on the other hand, we have 
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= [{s'^P^' - s^''p")p^'■' ■■■p^'- +■■■] (5) 
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with 



(for gi starting with n = 0, while for g2 we begin 
with n = 2). It has been argued that the n = 
moment of 32 should be [3] but it is not possi- 
ble to check this sum rule on a lattice. The lowest 
moment of gi is interesting because it can be re- 
lated to the fraction of the spin carried by the 
quarks in the nucleon, while g^ contains not only 
a„ (the so-called Wandzura-Wilczek contribution 
to 32 [4]), but also d„ - a twist-3 contribution. 
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Figure 1. At (vertical) -x (horizontal) section of 
the 4-dimensional lattice fitted onto the machine 
geometry. Each row of points represents one pro- 
cessor corresponding to one t value. Links are 
shown as lines. Every processor works simultane- 
ously on one link (e.g. the bold line). Even/odd 
lattice points are denoted by filled/shaded points. 



2. THE LATTICE CALCULATION 

The lattice geometry 

Wc perform our quenched QCD calculations for 
Wilson fermions on a N^NyN^Nt = 16^ x 32 lat- 
tice. For the gauge update we used a cycle con- 
sisting of a single 3-hit Metropohs sweep followed 
by 16 over-relaxation sweeps using the SU{3) al- 
gorithm as suggested in [5]. We repeated this 
cycle 20 times in order to generate a new con- 
figuration. The resulting configurations are high- 
ly independent - we see no correlations between 
hadronic quantities calculated on different gauge 
configurations. 

Our calculations are carried out on the 
Quadrics (or APE) Q16 machine, a parallel com- 
puter consisting of 128 processors arranged in 
a 2 X 2 X 32 periodic array. Since our lattice 
is 16^ X 32 each node stores a lattice region 
16 X 8 X 8 X 1. The fact that the block which 
each processor handles is only 1 lattice unit thick 
in the time direction means that considerable care 
is needed in programming the gauge update and 
fermion inversion. The Quadrics computer is a 
SIMD (Single Instruction, Multiple Data) ma- 
chine which means that every processor must per- 
form the same actions on its piece of the lattice as 
all the other processors. If the lattice were divid- 
ed up in the simplest possible way, as illustrated 
in Fig. 1, then there would be a conflict when 



updating space-hke links, since different proces- 
sors would be making simultaneous updates with- 
in the same plaquette, which would lead to false 
results. Our solution to this problem is illustrated 
in Fig. 2 in which we have used a new slanted co- 
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Figure 2. As for the previous Figure, but using 
slanted coordinates. 

ordinate system for dividing the lattice between 
the different processors. Thus we set 

Sx = {x -\- 1) mod 16 

sy = y 

Sz = Z 

St = t (7) 

and use in the programme always the slanted co- 
ordinates {sx, Sy, Sz, St). Now there are no con- 
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fiicts during the gauge update. An additional ad- 
vantage of the slanted coordinate system is that 
it allows algorithms that involve decomposing the 
lattice into even and odd sub-lattices. (Using the 
coordinate system of Fig. 1 this would not work - 
half the processors would be working on even sites 
and half on odd.) The even-odd decomposition is 
important when we invert the fermion matrix to 
find hadronic propagators, without this decom- 
position the inversion would be slower. 

Correlation functions 

The method to calculate matrix elements as in 
eqs. (2), (5) is well known, see [6] (we follow [7]). 
We first compute 2-point and 3-point correlation 
functions defined by 

Cv{t,T-p,0) 

= J2^p,a{Ba{t;p}Oir)Bp{0;p)), (8) 

where O is the Euclideanised version of the op- 
erators in eqs. (3), (6). The derivative operator 

-> ■(- 
D =D — D IS defined by 
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Using transfer matrix methods we can show that 
the ratio of 3- to 2- point correlation functions 
can be written as 

Rit,T;pSXO) 



(10) 



with 

F{T,J) = hi.rlTNJN]. 

N = 74 - ip- j/Ep + rn/Ep 
and J defined by 
{p, s\0\p, s) = u{p, s)Ju{p, s). 



(11) 
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When calculating 3-point functions it is partic- 
ularly important that the baryon operator B 



should be carefully chosen to have a low over- 
lap with excited baryon states, in order to make 
the plateau region in r as broad as possible. As 
our basic baryon operator we use the 'C75' wave- 
function (with C = 7472 in our representation) 
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with two important improvements. First we use 
a variation of 'Wuppertal smearing' [8] namely 
'Jacobi smearing' [9] in order to have an extended 
proton operator. Thus each quark operator in 
eq. (8) is replaced by 



(14) 



n=0 



We found a suitable value of {Ns,Ks) to be 
(50, 0.21) which for our present (/?, k) value of 
(6.0,0.155) gave a rms radius of about 4, corre- 
sponding roughly to ~ 0.5fm - half the nucleon 
radius. Secondly we replace each spinor by: 
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|(1 + 74)V^, 
V'Kl+74). 



(15) 



We call the resulting projected wave-function the 
'non-relativistic' wave function [10]. This re- 
placement leaves quantum numbers unchanged 
but we would expect it to improve overlap with 
those baryons which have slow-moving valence 
quarks. Practically this means that for each bary- 
on Green's function we invert on a smeared local 
source and consider only the first two Dirac com- 
ponents (i.e. the 2x2 sub-matrix of the upper 
left hand corner of the full Green's function for 
our representation of the 7-matrices) . So we only 
have 2x3 inversions to perform rather than the 
usual 4x3 inversions. This saves considerable 
computer time in the (minimal residue) inversion. 
In Fig. 3 we compare several baryon propagators. 
In each step we see an improved overlap with the 
proton. The projection operator is particularly 
effective at reducing the amplitude of the unwant- 
ed backward propagating state. Plots of the ef- 
fective mass \n{C{t)/C{t+ 1)), Fig. 4, show good 
plateaus for the two smeared operators with the 
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Figure 3. Comparison between point C75 wave 
function (dotted line, using a smaller sample size 
of 100), smeared wave function (dashed line) and 
smeared, non-relativistic wave function (full line) 
baryon correlation functions. For clarity we have 
dropped the (small) error bars. 

expected proton mass 0.66(1) [11] (the wall source 
result given there would seem to be a little low). 
In particular we see that after a distance of about 
4 time units we have very little trace of an excit- 
ed state. In Fig. 5 we show the nucleon energy 
(0.77(2)) with momentum p = pi = (27r/16, 0, 0). 
We see that the continuum energy dispersion re- 
lation = + m? is obeyed. 

To calculate 3-point functions we require addi- 
tional Green's functions, one for each chosen t, 
p and F. We have fixed t at 13, set p = 0, pi 
and have chosen two F matrices, F = ^(1 -I- 74), 
corresponding to the unpolariscd case, and F = 
1(1-1- 74)47572, corresponding to polarisation (-1- 

) in the y direction. We have also considered 

'tjj = u, d separately. This means that we must 
find 2 X 2 X 2 = 8 (half) Green's functions. We 
should emphasise here that the inserted operators 
O are neither smeared nor is a non-relativistic 
projection used (which means that the Green's 
functions running to O are smeared/projected 
on one Dirac index only). The choice f = 13 
is sufficient, larger values of t lead to unaccept- 
ably large errors in the signal for R - particularly 
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Figure 4. Effective mass for the point C75 wave 
function (empty square symbols) and the smeared 
non-relativistic wave function (filled circles). 
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Figure 5. Effective energy for the smeared non- 
relativistic wave function with momentum. 



for non-zero momentum. For example, test runs 
for t = 17 turned out to have ^ (9(2) larger er- 
rors, which roughly corresponded to the increase 
in the noise in the baryon correlation function 
from t = 13tot = 17. 

We have computed all operators up to and in- 
cluding 3 derivatives. At present we have generat- 
ed 350 configurations. We hope for a much higher 
statistic and so the following results should be re- 
garded as preliminary. 
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Mixing, renormalisation 

We must now choose favourable combinations 
of (T, /i indices for the operators in eqs. (3), (6). 
Simphcity dictates that wc choose indices where 
there is no mixing between the operator and low- 
er dimensional operators. On the lattice this is 
a more severe problem than in the continuum. 
We have explicitly checked [12] by seeing how the 
operator transforms on the lattice, whether mix- 
ing occurs or not. For low derivative operators 
there is no problem while for the higher deriva- 
tive operators mixing is excluded for certain off- 
diagonal combinations of the indices. However 
this then requires the use of non-zero momentum 
[13]. Note also that for the quenched theory there 
is no mixing with gluon operators [14]. 

The bare lattice operators are related to finite 
operators renormahsed at the scale /i by 



Oiti)=Zo{aii,gia))0{a). 
We define 



(16) 



{q{p)\0{nMp)) = {q{p)\0{aMp)) \%=,^, (17) 

where \q{p)) denotes a quark state of momentum 
p, and similarly for the gluonic operators. In 
the limit a — > this definition amounts to the 
continuum, momentum subtraction renormalisa- 
tion scheme. The structure functions themselves 
do not depend on //: the yu-dependence of the 
renormalisation constants is compensated by the 
//-dependence of the Wilson coefScients. If the 
Wilson coefficients c{Q'^ / fj? , g{ii)) and the struc- 
ture functions are evaluated at the momentum 
scale 



-2 2 

a = 1^1 



(18) 



(« 2GeV'^ here) the a- and //-dependence is elim- 
inated at least on the perturbative level. In the 
following we will assume (18), and we will denote 
the 'finite' part of the renormalisation constants 
by Zo without argument. 

For a few cases relevant to our investigation the 
renormalisation constants have been computed in 
perturbation theory. What is known is the axial 
current oq, where we shall take = 0.87 [15] 
and (a;), Z/^^^^ = 1.0127 [16]. We are in the pro- 
cess of calculating all Zo, see sect. 4. For the 



present in all other cases we take Zo = 1- (We 
expect that for 0, 1, > 1 derivatives Zo < 1, « 1, 
> 1 respectively [13].) 

3. RESULTS 

Wc first consider the unpolariscd structure 
functions. As noted previously, moments are di- 
rectly related to (a;"). For the lowest moment, 
(a;), we have no mixing problems, whether we 
choose off-diagonal components, 0{i4}, method 
(a), or diagonal components, 044 — (On -I- O22 + 
C33)/3, method (6), which provides another use- 
ful consistency check (and perhaps gives indica- 
tions about possible lattice artifacts). For (a) we 
have from cq. (10), i?a = i/^(a;)„-l/2K-pi-(a;)a. In 
Fig. 6 we show Ra for the u quark. From eq. (10) 
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Figure 6. The ratio of the 3-point connected cor- 
relation function to the 2-point correlation func- 
tion, R^a\ versus T for {x)a'\ A linear fit, dotted 
line, m 4 < r < 9 gives 0.53(5). For R^a^ we have 
0.23(2). 

we hope that in the region ^ r ^ t = 13 wc sec 
a plateau. In addition the dotted line is the 'con- 
jugate signal' (here the real part of the correlation 
function). It should be zero, and we monitor it 
to give a further indication of how good our sig- 
nal is. We would claim that we see a reasonable 
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signal. A linear fit gives 

« 0.42(4), « 0.18(2). (19) 

Similarly for (6) wc have R = —l/Zi^^^^ ■ 1/2k ■ m ■ 
(x) shown in Fig. 7, with results 
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Figure 7. R^^^ versus t for A linear fit 

gives -0.99(3). For R^f' we have -0.45(3). 

(a:)!"' « 0.47(2), (x)*^' « 0.21(2). (20) 

The methods seem to be consistent. The numbers 
can be compared to plionomenological parame- 
tcrisations of the experimental data, e.g. from 
[17] (L»_), of 0.28, 0.11 at Q'^ = iGeV'^ for u, 
d respectively. Analytically continuing in ra, for 
{x^) we have with O^n^ - (£'{224} + C'{334})/2, 
R = -l/Z(^^2) ■1/2k-p\- (x^) shown in Fig. 8, 
which gives 

(a;2)(") = 0.12(2), {x'^)^'^') = 0.050(8) (21) 

to be compared with 0.08, 0.03 respectively. Both 
the (a;), (a;^) numerical results seem to be about 
50% larger than the phenomenological values as 
already found for larger quark masses in [6]. Fi- 
nally we have computed (x^). With O^imy + 

C{2233} - C'{1133} - C'{2244}, R = 1/^(.t3> ' 1/2k • 

^PiPi ■ (^^) shown in Fig. 9. Although there 
are large fluctuations in the data, the values 

= 0.031(8), (a;^^)^'*) = 0.0076(47) (22) 



0.08. 
0.06. 



0.04. 
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Figure 8. i?^^^ versus r for {x'^)^^\ A linear fit 
gives -0.062(8). ForR^'^^ we have -0.025(4). 



are consistent with the phenomenological results 
of 0.03, 0.009. 

We now turn to the polarised cases and start 
with flo , which is the matrix element of the axial 
current. In Fig. 10 we show i?^"' for 0^-2 which 
gives i?, = i/Zag ■ 1/2k ■ rn/2Eff ■ uq with p = 0. 
Conventionally we set oq = 2Aq. A Hnear fit gives 

Ag("' « 0.86(5), Ag^"*) « -0.27(2). (23) 

These numbers are to be compared with the 
experimental results, e.g. as given in [18] of 
0.80(4), -0.46(4) respectively and -0.13(4) for 
the strange quark. For the d quark, at least, 
the result seems a little large. Non-singlet quan- 
tities are likely to be more independent of sea 
quark contributions (which are absent in the 
quenched approximation). For the axial decay 
constant (from /3 decay) we have experimentally 
Am — Ad = K 1.26. Our result « 1.13 seems 
a little small compared to this. 

For comparison we show in Fig. 11a calculation 
of the same quantity, but with non-zero momen- 
tum. Wc see a degradation of the signal by a 
factor of roughly 0(2). (To obtain a comparable 
signal as for the zero momentum case would thus 
require 0(4) times more configurations.) Never- 
theless a consistent result with p = is obtained. 
We find 0.88(12), -0.31(6) for u, d respectively. 



7 




Next we consider 02 . To avoid mixing problems 
we choose the operator 05;{2i4} so the matrix el- 
ement needs p = pi- We have R = i/Za^ ■ 1/2k ■ 
1/6 ■ mpi ■ 02 shown in Fig. 12. As expected due 
to the complexity of the operator, the signal is at 
present not so good. Fits lead to a'^^ « 0.14(3), 



-0.0047(12), which gives 



0.016(3) 
0.0034(10) 



(proton) 
(neutron) 



(24) 



(for the neutron u ^ d). This is to be compared 
to the EMC result of 0.022 for the proton. 

For d2 we choose 05;[2{i]4} and show the ratio 
R = ■ 1/2k • 1/3 • mpi ■ d2 in Fig. 13. An 

estimate gives ^2 ~ —0.18(2) for the u quark and 
« 0.0036(11) for the d quark. Thus from eq. (4) 
we have 



j dxx^g2 « I 



-0.013(2) 
-0.003(1) 



0.010(2) 
0.002(1) 



(25) 



(proton/neutron) where the first number comes 
from 02 (twist 2), while the second comes from ^2 
(twist 3). They arc of the same order of magni- 
tude. Two other theoretical estimates have been 
made, using bag models [19] and sum rules [20]. 
Sum rule results suggest that for the proton ^2 is 
very small. At present we seem to have a some- 
what larger result. 



4. RENORMALISATION 

Parallel to our numerical work we have begun 
to compute the renormahsation constants for ah 
twist-2 operators up to spin 4, to one loop order. 
Though this is a weh defined task, it is compu- 
tationally voluminous due to complications im- 
posed by the finite periodic volume and the hy- 
percubic symmetry. In particular the vertices as- 
sociated with the higher dimensional operators 
are of unprecedented complexity. To master this 
task we are currently developing packages of com- 
puter algebraic programmes using Mathematica 
and Maple. 

Gluon Operators 

The operator we have considered first is the 
gluon operator 



(ng — TrF F 



where 



" tJl,,v—± 



(26) 



(27) 



and Ufj,!, is the plaquette operator. This operator 
gives us the lowest moment of the gluon distribu- 
tion. We need to expand F^,^ in terms of powers 
of the coupling constant. Writing [21] 

U^, = exp[j(ff,^W + + ff3,^(fj + ■ ■ ■)], (28) 
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Figure 11. i?'"' versus r for a''^'^ with momen- 
tum. A linear fit gives 2.8(4). For R^'^^ we have 
-0.98(19). 

where (p^S is of order k in the gauge potential, we 
obtain 



1 

4^ 



+ 



(29) 



The calculation of the loop integrals is done in 
two parts. Following [22] we write 

Hp) {IhihAp) - IhihAp)) +IniH2{p), (30) 

where 



hinAp) 



Pai ' ' ' Pan 

n=0 



! ■ ■ ■ 



-IniH2{lJ)\p^o- (31) 



The first term in eq. (30) is ultraviolet finite and 
is computed in the continuum, while the second 
term is ultraviolet divergent and is computed on 
the lattice. Both terms have infrared divergences. 
They arc regulated by dimensional rcgularisation 
and cancel out in the sum. The calculations are 
done in Feynman gauge. The programme has 
been developed to such a level that all what is 
needed as input is to state the Feynman rules in 
symbolic form, both for the continuum and the 




6.000 
/t=0.1550 



Figure 12. versus r for 02"^ . A linear fit 

gives 0.020(4). For R'^^^ we have -0.0066(17). 

lattice part of the calculation. The lattice inte- 
grals are reduced to a few integrals that can be 
computed to arbitrary precision. If we write 



l-g^N,Bg, 



find 



B„ = - — 



3_ 
64 



487r2 



^ Z ^ Z + ^ 



(32) 



(33) 



where Zq = 0.15493339, Z^ = 0.10778131. This 
is, as all our results, in the quenched approxima- 
tion and agrees with [16]. At f3 = 6.0 this gives 
Zg = 1.296. The calculation can be extended 
without any difficulties to higher moments of the 
gluon distribution function. 

Quark Operators 

The calculation of the renormalisation con- 
stants for the quark operators is in progress. Here 
the main problem is to compute the lattice inte- 
grals over the fermion propagators to a sufficient 
accuracy. Like in the gluonic case we found that 
it is possible to reduce all integrals to just a few 
basic ones. We hope to report our results shortly. 

5. DISCUSSION 

We have presented here an outline of our pro- 
gramme for calculating low moments of polarised 
and unpolarised structure functions. At present 
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Figure 13. versus r for ■ A linear fit 

gives -0.049(6). For R'^'^^ we have 0.010(3). 

we have only one k value, which represents rather 
a heavy quark mass of about 66MeV. We need 
at least one more k value (two would be bet- 
ter) to be able to extrapolate to a more phys- 
ical quark mass. Finally, although also not re- 
ported here, we are monitoring the gluon mo- 
ments (from FF and FDDF) and calculating 
the hadron mass spectrum, electro-magnetic form 
factors (vector current), the scalar (sigma term) 
and pseudoscalar density (for the Qnn-k coupling), 
all at non-zero momentum transfer [12]. 
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